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SUMMARY 7975 62^- 


We develop a new numerical approach to study the spatially evolving instability of the 
/ streamwise dominant flow in the presence of roughness elements. The difficulty in handling the flow 
! over the boundary surface with general geometry is removed by using a new conservative form of 
the governing equations and an analytical mapping. The numerical scheme uses second-order 
backward Euler in time, fourth-order central differences in all three spatial directions, and 
boundary-fitted staggered grids. A three-dimensional channel with multiple two-dimensional-type 
roughness elements is employed as the test case. Fourier analysis is used to decompose different 
Fourier modes of the disturbance. The results show that surface roughness leads to transition at 
lower Reynolds number than for smooth channels. 


INTRODUCTION 


Transition from laminar to turbulent flow is a phenomenon of great importance and practical 
interest. Major experimental work in aerodynamics has been pursued to study boundary layer 
transition, and this in turn has led to a critical need for further understanding of this fundamental 
process. Unfortunately, to date, no reliable methods exist for predicting transition in the presence of 
surface roughness, either experimentally or numerically. In fact, there are many factors that affect 
transition, such as solid wall temperature, solid wall curvature, pressure gradients, free-stream 
disturbance, and surface roughness. There has been some experimental activity dealing with the 
effect of both 2-D and 3-D roughness elements on transition from laminar to turbulent flow. For 
example, Klebanoff et al showed that surface roughness induces early transition [1]. Also, others 

have obtained limited numerical results with 2-D flow [2]. 

The purpose of the present work is to develop new efficient and easy-to-use methods for 
numerical simulation of the effect of surface roughness on flow transition. A new conservative form 
of the governing equations is derived. Because of the high sensitivity of transitional flows, a high 
order scheme based on our earlier work (cf. [3] - [6]) is developed. For the grid generation scheme, 
an analytical map is used, so the Jacobian coefficients are computed exactly. Moreover, we develop 
the governing equation in a form that enables a much simpler numerical process. 

We impose single and multiple 2-D-type roughness elements on the lower solid wall to test the 
effect of surface roughness on flow transition. A Fourier transformation is employed to analyze 
different modes of the resulting disturbance. The computational results show that the induced 
mean flow distortion and other high frequency waves make the flow more unstable. 

* This work was supported by NASA under grant number NAS1-19312. 
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GOVERNING EQUATIONS 


In this study, the three-dimensional, time-dependent, incompressible Navier-Stokes equations, 
which axe nondimensionalized by the channel half-height h and the centerline velocity U . are 
considered as the governing equations for 3-D channel flow (Figure 1): 
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where u, v, and w are velocity components in the x-, y- } and z- directions, respectively, P is the 
pressure, and Re is the Reynolds number based on the centerline velocity U K of mean flow, the 
channel half-height h, and the viscosity parameter v. 


v 



Figure 1. 3-D channel with a single 2-D-type roughness element. 
For the current work, we consider a special mapping 
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which implies that 



and the final forms of the momentum and continuity equations are: 
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for the perturbation flow, where, 


A > = ^5 + ( nl + ’ll + + ^5 + 2 + 2».^ + <x~ + a»» + 

The inverse transformation of the variables under this special mapping becomes 
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Here we have seven unknowns (u, v, w, P, U> V, W), seven equations ((6) - (9), (15) - (17)) for 
the base flow or total flow, and seven equations ((10)-(13), (15) - (17)) for the perturbation. 

Our solution process is outlined as follows: 

1. Perform the surface and grid generation process to obtain the required Jacobian coefficients. 
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2. Solve system (6) - (9) and (15) - (17) to obtain the base flow solution. 

3. Solve system (10)— (13) and (15) - (17) to obtain the perturbation solution based on the above 
base flow. 


For the channel flow, though the boundary conditions axe quite simple, there are still some 
difficulties we need to overcome. The so-called buffer domain [7] technique is used here for both the 
base flow and perturbation. For details, see [6]. The boundary condition for the solid wall is the 
no-slip boundary condition: 

Inwall vail ^ wall ^ wall Vt vail ^tuoiJ 0* (l A) 


No boundary condition is needed for the pressure at the solid wall since we use a staggered grid. 
For the inflow, Poiseuille flow is imposed at the inlet for the base flow solution, and the 
eigenfunctions obtained from the linear stability theory with specified Reynolds number axe 
employed at the flow inlet for the perturbation. The final outflow boundary conditions are: 
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• for the perturbation flow, 
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and the associated u , v, w can similarly be obtained by the inverse transformation (15) - (17). 


Periodicity is assumed in the spanwise direction. 


NUMERICAL PROCEDURE 
Surface and Grid Generation 


Assume that no stagnation points exist in the computational domain. Solitary type roughness 
elements are overlapped on the lower solid wall of the channel to simulate surface roughness. For 
the 2-D roughness elements (because the grid is uniform and the domain is periodic in the spanwise 
z direction, we need only discuss the 2-D case here), the surface can be expressed as 

m 

f( x ) = Y. Kisech 2 {b t (x - x,)), (21) 

i=i 


380 


where Ki is the height of the roughness element, 6; is a parameter for adjusting the curvature rate of 
the roughness element, and Xj is the peak point coordinate. 



Figure 2. Physical, stretched, and uniform grids. 

Our grid generation approach consists of two mappings(see Figure 2): 

• from the physical grid that conforms to the rough boundary to the stretched intermediate grid 
that is uniform in the x direction but nonuniform in y, 


• from the stretched grid to the uniform computational grid. 


The resulting mapping from the physical to the uniform computational grid is 
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where y max is the maximum height of the computational domain and a is the parameter for 
adjusting the density of grids near the lower solid wall. The required Jacobian coefficients sure 
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Discretization 


In the computational (£, 77, () space, the grids are uniform. Suppose u, v, w and U , V, W are 
defined in terms of a staggered grid in the computational space (see Figure 3). Here, the values of P 
are associated with its cell centers, u and U with centers of the cell surfaces parallel to the (77, () 
plane, v and V with centers of the cell surfaces parallel to the (£, £) plane, and w and W with 
centers of the cell surfaces parallel to the (£, 77) plane. 

Second-order backward Euler differences are used in the time direction, and fourth-order central 
differences are used in space. Details of such an approach can be found in [5]. With those 
assumptions, we can write the discretized governing equations symbolically as follows: 

Aee u ee + Aeue + Awuw + Aww^ww + Annum* + Anun + Asu$ + Assess + 


Aff^ff + Apup + Agujg + Abb^bb — Acuc + DwwPww + DwPw + DePe " DcPc = S u , ( 29 ) 
Beevee 4- Beve + Bwvw + Bwwvww + Bnnvnn + Bnvn + -Ss v s + Bssvss + 

Bff v ff + Bpvp + BgVB + Bbb v bb — Be vc + EssPss + EsPs + EnPn — EcPc = S v , ( 30 ) 
Cee w ee + Cevje + Cwww + Cwwvjww + Cnnvjnn + Cnwn + Csvis + Csswss + 

Cffwff + Cfwf + Cbvib + Cbbvjbb — Ccvjc + FbbPbb + FbPb + FfPf — FcPc = S w , ( 31 ) 
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DV S V S - DV c Vc + DWffWff + DW f W f + DW b W b - DW c Wc = S m , ( 32 ) 

u c = Vy C U c , ( 33 ) 

v>c = Vy C W c , ( 34 ) 

vc = v: c U c + Vc + v: c W c . (35) 



Figure 3. Staggered grid structure in the computational (£, 77, Q space. 


The coefficients and source term for the interior points of the discrete £— momentum equation (29) 
associated with uc are given as follows: 
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Here, superscripts n and n — 1 are used to indicate values at previous time steps, and the 
superscript n 4- 1, which indicates the current time step, is dropped for convenience. Lower case 
subscripts denote the approximate values of the v and w at points where the associated values of u 
are located (Figure 4). Other symbols used in the above formulas are as follows: 
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Figure 4. Neighbor points for £— momentum equation 
(U are at the same points as u and are not shown here). 


All function values that are required at other than the canonical locations are obtained by 
fourth-order interpolation in the computational space. For example (see Figure 5), 

V c — (9 (Vc + v N + Vnw + Vw) — (Vsww + Vnnww + Vse + Vnne))/ 32, (37) 

P a — (9 (Ps + Psw) ~~ (PsE + PsWw))/l-6 • (^) 

The coefficients for the t ) — and £ — momentum equations are defined in an analogous way, the 
discrete continuity equation is developed simply by applying central differences to each terms. 

On the solid wall boundary points, we change the rj— direction difference to second order, and 
maintain fourth-order in both the £— and directions. For more details, see [6]. 
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Figure 5. Neighbor points for fourth-order approximation for V c and P a . 
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Line Distributive Relaxation 


We use here the same basic line distributive relaxation method developed in our previous work 
(cf. [5]), with some modifications. Figure 6 shows the distribution of corrections for the group of 
variables that axe located in the (£,77) plane. 

The process that we use for solving the discrete system (29)-(35) can be described as follows: 

• Freezing P, U, V, W, v, and w, perform point Gauss-Seidel relaxation on (29) over the entire 
computational domain to obtain a new u. 

• Freezing P, U, V, W, it, and w, perform point Gauss-Seidel relaxation on (30) over the entire 
computational domain to obtain a new v. 

• Freezing P, U, V , W, it, and v, perform point Gauss-Seidel relaxation on (31) over the entire 
computational domain to obtain a new w. 

• Use transformation (33)-(35) to obtain new U, V , W. 

• For all j = 2, 3, * ■ * , Tij 1 at once: change U%jky Ut+i jky Vijky fc+i fo satisfy - the 

continuity equations, then update P^* so that the new U, V, W and P as well as the 
associated transferred u, v, w satisfy the three momentum equations. 


£ 

Figure 6. Distribution of corrections in the (£, rf) plane. 
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Remark. Since all of the u, v, w have been previously relaxed, and the U, V, W are updated to 
perform the latter corrections, we assume that equations (29)-(31) hold exactly. Let e, 6, a and 7 
represent the corrections for U , V, W and P , respectively. Thus, for cube ijk, when we distribute 
the corrections according to Figure 6, the correction equations corresponding to (29)-(31) are 

(Afrtf* • + A% k W )cj - Df-y, = 0, (39) 

- <5, +1 ) - - <■>,-) - = 0. (40) 

(CgV’**’ + Cfi/?*)*) ~ = 0. (41) 

(Duf + Dirpy, + (Dwp + Dwfy 7, + Dv^\s j - s j+l ) - nvg*^., - s t ) = sp, (42) 

j = 2 , 3 ,- •• ,n 5 - 1. 

This system has 4 (rij — 2) equations and 4 (rij — 2) variables. Unfortunately, coupling between the 
correction variables makes the problem somewhat complicated. To develop a simpler approximate 
system, define 
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With the above approximations, uj xj and u) ZJ can be treated as known parameters, so equation (44) 
can be written in terms of the unknowns 6j only: 
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Then we obtain the tridiagonal system 
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Thus, 6j, j = 2, 3, • • • , rij — 1 can be determined very efficiently. The other velocity corrections are 
given by 


Cj — ui x j6j, <7j — u z j6j 1 

3 = 2j 3> • ■ • , Tij 1. 

The u, v, and w are then updated on all cells in the i, k y— line as follows: 
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Finally, the pressure corrections 7 j are determined as follows: 
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COMPUTATIONAL RESULTS 


We first tested the code by applying it to single and double roughness elements using a moderate 
sized grid. A 170 x 50 x 4 grid, including a five T-S wavelength physical domain and a one 
wavelength buffer domain, was used. Considering that Re = 5000 corresponds to a decreasing mode 
according to the linear stability theory, we use this Reynolds number in our code to investigate the 
effect of roughness elements. Fourier analysis is used to decompose different Fourier modes of the 
disturbance. For details about the Fourier transformation approach, see [8]. 

Let kj = 0.15 and b t = y/2. For the single roughness case, the peak point is at i = 30; for the 
double roughness case, the peak points are at i = 42 and i = 70. The stretch parameter a is set to 
4, and the amplitude of the disturbance e is set to 0.0025v 2. Figure 7 displays contours of the 
perturbation streamfunctions, showing that the roughness elements make the disturbance increase 
for a certain distance downstream. The results of Fourier transformation given in Figure 8 show 
that the mean flow distortion and first and second harmonic waves are amplified over this distance. 

To test the effect of multiple elements, a 402 x 66 x 4 grid (including a nine wavelength physical 
domain and a one wavelength buffer domain) is used. Here we set Kj =0.12, bi = 2.0, and a = 4.5. 
The first roughness is at i = 82, which is two wavelengths from the inflow boundary. We placed 
seven roughness elements in the computational domain, starting from i — 82 and spared 40 grid 
points apart. Figures 9 and 10 depict the contour plots of streamfunctions and vorticities, showing 
very clearly that the disturbance is amplified after it passes each element. 




(b) 


Figure 7. Contours of perturbation streamfunctions with Re = 5000, /t; = 1.5, 
e = 0.0025\/2 and 170 x 50 x 4 grid. Flow direction is from left to right. 
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10 z rms 










Figure 8. Maximum amplitudes of fundamental wave Ui, Vi, mean-flow distortion u 0 , v Q , 
first harmonic wave it2, U2, and second harmonic wave U3, C3 for Re = 5000, 

Ki = 0.15, and e = 0.0025\/2 with two roughness elements (grid: 170 x 50 x 4). 
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(b) 


Figure 9. Contour plots of perturbation streamfunctions and vorticities 
for Re — 5000, Ki = 0.12, and e — 0.0025v^ with seven roughness 
elements (grid: 402 x 66 x 4, flow direction is from left to right). 


total atreamfunction contours 



(a) 

spawnwise total vorticity contours 



(b) 


Figure 10. Contour plots of total streamfunctions and vorticities for 
Re = 5000, Ki = 0.12, and e = 0.0025\/2 with seven roughness 
elements (grid: 402 x 66 x 4, flow direction is from left to right). 


CONCLUDING REMARKS 


As expected on physical grounds, we find that the spatial growth rates of the disturbance 
increase when surface roughness is present. Though our work is limited to roughness without 
stagnation points in the computational domain, such a scope includes a rather large variety of real 
roughness elements. Moreover, the code is very efficient, requiring about 2.68 seconds per time step 
for the 402 x 66 x 4 grid case (equivalent to about 26 [is per time step per grid point). 
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